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Abstract 

The /-plane and /3-plane wave propagation properties are examined for discretisations of the 
linearized rotating shallow- water equations using the P\r,Q-P2 finite element pair on arbitrary 
triangulations in planar geometry. A discrete Helmholtz decomposition of the functions in the 
velocity space based on potentials taken from the pressure space is used to provide a complete 
description of the numerical wave propagation for the discretised equations. In the /-plane (planar 
geometry, Coriolis force independent of space) case, this decomposition is used to obtain decoupled 
equations for the geostrophic modes, the inertia-gravity modes, and the inertial oscillations. As 
has been noticed previously, the geostrophic modes are steady. The Helmholtz decomposition is 
used to show that the resulting inertia-gravity wave equation is third-order accurate in space. In 
general the P1dg-P2 finite element pair is second-order accurate, so this leads to very accurate 
wave propagation. It is further shown that the only spurious modes supported by this discretisation 
are spurious inertial oscillations which have frequency /, and which do not propagate. A restriction 
of the PIdg velocity space is proposed in which these modes are not present, leading to a finite 
element discretisation which is completely free of spurious modes. The Helmholtz decomposition 
also allows a simple derivation of the quasi-geostrophic limit of the discretised P\dg-P^ equations 
in the /3-plane (planar geometry, Coriolis force linear in space) case resulting in a Rossby wave 
equation which is also third-order accurate. This means that the dispersion relation for the wave 
propagation is very accurate; an illustration of this is provided by a numerical dispersion analysis 
in the case of a triangulation consisting of equilateral triangles. 

Keywords: Mixed finite elements, geophysical fluid dynamics, Rossby waves, spurious modes, 
numerical weather prediction 
2010 MSC: 65M60 



1. Introduction 

Recently there has been growing interest in developing more general horizontal discretisation 
schemes for numerical weather prediction (NWP) models with computational meshes constructed 
from triangles or hexagons. There are two principal motivations for this. Firstly, geodesic grids 
(which are obtained by iterative refinement of an icosahedron using triangles, sometimes trans- 
forming to the dual grid which is a mesh of hexagons with exactly 12 pentagons located at the 
vertices of the original icosahedron) provide similar grid cell areas over the entire sphere, which has 
possible advantages for accurate representation of wave propagation. Furthermore, geodesic grids 
also avoid the very fine grid cells obtained near the North and South poles on latitude-longitude 
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grids, which lead to large Courant numbers, and cause bottlenecks in communication between 
processors on parallel systems. This has led to a number of groups developing weather and climate 



models which use geodesic grids (Ringler et al. , 2000 Majewski et al. , 2002 Satoh et al. 2008). 



Secondly, triangles facilitate the implementation of adaptive mesh refinement. This allows nested 
regional models within a global model, and further allows dynamic mesh refinement in which the 
mesh resolution is locally modified in response to the dynamics in the course of a forecast. The 
development of new numerical schemes that correctly represent the qualitative properties of wave 
propagation on these grids, and under adaptive-mesh refinement, is crucial. 

Possible discretisations on triangular or hexagonal meshes are obtained using three different 
approaches: finite difference methods, finite volume methods and finite element methods. To 
eliminate spurious pressure modes, finite difference methods use a C-grid in which the edge-normal 
velocity is stored at the edge-centres, and the pressure is stored at the cell-centres. On quadrilateral 
grids, the wave propagation is observed to be well represented provided that the Rossby radius is 
well- resolved (Arakawa and Lamb 1977 Fox-Rabinovitz 1996; Randall 1994). On triangular and 
hexagonal grids the problem lies in finding a scheme for reconstructing the Coriolis force (which 
requires the tangential velocity) from the normal velocity. Recently, a reconstruction scheme 
was found which results in steady geostrophic modes for C-grid discretisations on the regular 
hexagonal grid in the plane (Thuburn, 2008). In the same paper it was shown that the resulting 
discrete system on the /3-plane has a spurious extra Rossby wave branch, with very slow Eastward 



phase velocities. This reconstruction was extended to arbitrarily structured C-grids in Thuburn 



et al. (2009). The finite element method provides a great degree of flexibility in the choices of 
discretisation for velocity and pressure. Amongst the many finite element pairs that have been 
proposed for the rotating shallow-water equations are the PItvc-PI and Pl-iso P2-P1 elements 
(investigated and compared to several other element pairs in Le Roux et al. (1998)), the RT0 



elements (introduced in Raviart and Thomas (1977) and proposed for the shallow- water equations 
in Walters and Casulli|( 1998 )) and equal-order elements with stabilisation (also proposed in Walters 



and Casulli (1998)). 



In this paper we study the wave propagation properties of the recently proposed PI dg~P2 finite 
element discretisation. This discretisation uses a mixed finite element pair with The PI dg-P2 finite 
element discretisation was introduced in Cotter et al. (2009b), and was designed to accomodate the 
geostrophic balance relation between pressure and velocity without introducing spurious pressure 
modes. This is achieved by using a quadratic (P2) continuous finite element basis for pressure, 
and a linear discontinuous (PIdg) finite element basis for velocity (hence the name P1dg-P2). 
The pressure polynomials are one order higher than the velocity polynomials, which accomodates 
the geostrophic balance relation since the pressure gradient and the velocity are both linear within 
each element. Making the velocity basis discontinuous increases the number of velocity degrees of 
freedom so that there are no spurious pressure modes. The lack of pressure modes was investigated 
numerically in Cotter et al. (2009b) and subsequently proved in Cotter et al. (2009a), where it was 
also shown that this combination of spaces means that geostrophically balanced states are exact 
steady states of the linear equations on arbitrary unstructured meshes (this property can also be 
obtained for C-grid finite difference methods as described in Thuburn et al. (2009), with the added 
restriction that the meshes satisfy an orthogonality property). In this paper we go further and 
produce a complete description of the numerical wave propagation properties of P1dc?-P2, which 
is facilitated by the construction of a discrete Helmholtz decomposition of the PIdg space. 

The rest of this paper is organised as follows. In Section [2j we show that P\j^q-P2 has a 
discrete Helmholtz decomposition. In Section [3] we use this decomposition to analyse the wave 
propagation on the /-plane. We show that there are three types of modes: steady geostrophic 
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modes, inertia-gravity modes, and inertial oscillations (of which only one is a physical mode). 
We show that the inertial oscillations do not propagate and can be filtered out by solving two 
discretised elliptic equations. We also show that the velocity may be eliminated to obtain a third- 
order accurate inertia-gravity wave equation, and hence claim that the wave propagation is very 
accurate on arbitrary unstructured meshes. In Section |4| we use the Helmholtz decomposition to 
analyse the Rossby wave propagation on the /3-plane in the quasi-geostrophic limit (following the 
approach of Thuburn (2008)). We obtain a third-order accurate Rossby wave equation, and hence 
claim that the Rossby wave equation is also very accurate. Finally, in section [5] we give a summary 
and outlook. 



2. Discrete Helmholtz decomposition for P\dg-P1 

In this section we show that the PI dg~P2 finite element discretisation has a discrete Helmholtz 
decomposition for PIdg-P 1 ^- We shall adopt the notation that the 5 superscript indicates a 
numerical approximation in a finite element space; functions without subscripts indicate continuous 
fields. We start by stating two properties of P1dg-P2 which we shall use throughout. 

Definition 1 (Embedding conditions). Let V be the chosen vector space of finite element velocity 
fields (in the case of P1dg-P2, V is the space PIdg of velocity fields u s that are linear in each 
triangular element, with no continuity constraints across element boundaries), and let H be the 
chosen vector space of finite element pressure fields (in the case of P1dg-P2, H is the space P2 of 
pressure fields h s that are quadratic in each triangular element and are constrained to be continuous 
across element boundaries). 

1. The operator V defined by the pointwise gradient 

q\x) = Vh s (x) 

maps from H into V. 

2. The skew operator _L defined by the pointwise formula 

q 5 (x) = (u s (x)) ± = (-uiu s 1 ) 

maps from V into itself. 

These are the only conditions that we use in the paper and hence any properties extend to any 
other finite element pair that satisfies these conditions (P0-P1 or PnoG-P{ n + 1) with any n > 1, 
for example). 

These conditions are most definitely not satisfied by all possible pairs (V, H), as illustrated by 
the following examples. 

Example 2 (Pl-Pl). The finite element pair known as Pl-Pl (which may be used for the shallow- 
water equations but requires stabilisation as described in Walters and Casulli (199®) ) is defined as 
follows: 

• The mesh Ai is composed of triangular elements. 

• H is the space of elementwise-linear functions h s which are continuous across element bound- 
aries. 

• V is the space of vector fields u 5 with both of the Cartesian components (u 5 ,v 5 ) in H. 
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Condition 1 of Definition^ is not satisfied by the Pl-Pl pair since gradients of functions in H 
are discontinuous across element boundaries. Condition 2 is satisfied since the same continuity 
conditions are required for normal and tangential components. 



Example 3 (RTO). The lowest order Raviart- Thomas (Raviart and Thomas, 1971) velocity space 
(known as RTO) is constructed on a mesh M. composed of triangular elements. It consists of 
elementwise constant vector fields which are constrained to have continuous normal components 
across element boundaries. RTO does not satisfy condition 2 of Definition^ since the _L operator 
transforms vector fields with discontinuities in the tangential component (which are permitted in 
RTO) into vector fields with discontinuities in the normal component (which are not). 

We now describe some examples of finite element pairs which do satisfy the conditions in 
Definition Q] 

Example 4 (P0-P1). The finite element pair known as P0-P1 (applied to ocean modelling in 
Umgiesser et al. (2004\ ), and analysed in Roux et al. (2001); Roux and Pouliot (200$fy ) is defined 



as follows: 

• The mesh Ai is composed of triangular elements. 

• H is the space of elementwise-linear functions h s which are continuous across element bound- 
aries. 

• V is the space of elementwise- constant vectors with discontinuities across element boundaries 
permitted. 



Example 5 (P1dg-P2). The finite element pair known as PIdg~P2 (Cotter et al, 2009b) is 
defined as follows: 

• The mesh M. is composed of triangular elements. 

• H is the space of elementwise- quadratic functions h s which are continuous across element 
boundaries. 

• V is the space of elementwise-linear vectors with discontinuities across element boundaries 
permitted. 

Each of these examples satisfy both conditions in Definition [TJ condition 1 holds because taking 
the gradient of a elementwise polynomial n — 1 which is continuous across element boundaries 
results in a vector field which is discontinuous across element boundaries and is composed of 
elementwise polynomials of one degree n, and condition 2 holds since the velocity space uses 
the same continuity constraints for normal and tangential components e.g. both components are 
allowed to be discontinuous. This defines a whole sequence of high-order Pn£,G~P{n + 1) element 
pairs. Similar elements can be constructed on quadrilateral elements. Since we only require these 
two conditions to prove our optimal balance property which holds on arbitrary meshes, we can 
also construct finite element spaces on mixed meshes composed of quadrilaterals and triangles, 
for example. It is also possible to use p-adaptivity in which different orders of polynomials are 
used in different elements, as long as the conditions are satisfied. To make the rest of the paper 
less abstract, we shall only discuss P1dg-P2, but all of the results are easily extended (with th 
appropriate orders of accuracy) to any element pair satisfying Definition [T] 
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Next we note that the gradient and skew-gradient any two pressure fields , ip in the pressure 
space P2 are orthogonal in the L2 inner product, 

(v^ 5 ,vV> = / v^-vVdy = o, 

Jn 

where Q is the solution domain which is either the sphere, or periodic boundary conditions. This 



was proved by direct computation in Cotter et al. (2009a). Hence, any velocity field u in PIdg 



can be written uniquely in an orthogonal decomposition 

u 5 = u 5 + V(j) 5 + V x ip s + u s , (1) 

where u s is independent of space, where 4> s and ip s are both in the space P2, which consists of P2 
functions with mean zero, i.e. 

<0 5 ,1>= / <p 5 dv = o, <^,1>= / iP 5 dV = 0, 
Jn Jn 

and where u s is orthogonal to the gradient or skew-gradient of any P2 function or, i.e. 

(u 5 ,Va 5 ) = (u^^a 6 ) =0. 



Furthermore, if any such ii satisfies 



(u 5 ,u 5 ) = 0, 



then it 5 = 0, since ii s is obtained from orthogonal completion. In general the dimension of the 
orthogonal subspace containing the vector fields of the form ii is non-zero, since there are more 
than twice as many degrees of freedom in the velocity space V as the pressure space H. The 
dimension of V is 6nj (where rif is the number of elements), and the dimension of F is n v + n e 
(where n v is the number of vertices and n e is the number of edges). For doubly periodic boundary 
conditions, Euler's polyhedral formula on the torus then gives dim(H) = n v + n e = 2n e — nj. For 
a triangulation, 2n e — 3n/ since each triangle has three edges which are each shared between two 
faces, so dim(if) = 2nj < 3n/ = dim(V)/2. Since 2dim(if) < dim(V) it is not possible to span 
V entirely with functions of the form V^ip 6 + <fr s G H, and so components of the form ii s 

will always be present. 

Equation ([T]) is identical to the Helmholtz decomposition for arbitrary continuous velocity fields 
in which any continuous velocity field u can be written as a constant plus a gradient of a potential 
plus the skew gradient of a streamfunction; the only difference in the discrete PI dg-P^ case is the 
extra component ii 5 . This extra component gives rise to spurious inertial oscillations in the P1dg~ 
P2 finite element discretisation applied to the rotating shallow-water equations. It is possible to 
describe a reduced velocity space, which we call H(P2), consisting of velocity fields which can be 
written as 

v 5 = v s + V(j) 5 + V ± V 5 , 

where v s is independent of space, where <p s and ip s are both in the space P2, i.e. we have removed 
the spurious velocity component. It is possible to project a PIdg velocity field u s into H(P2), by 
first computing the mean component, 



u 



In^V ' 
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and then extracting the velocity potential and streamfunction by solving 

(Va s ,V(j) S ) = (Va 5 V), 

and 

<Vc/,W> 5 ) = (VWV), 

for all P2 functions a s . This amounts to solving elliptic problems for (p 5 and ip s . Then, the 
projection of u 5 into H(P2) is given by 

u s + V(j) s + v V 5 - 

3. Discrete wave propagation on the /-plane 

In this section we describe all of the numerical solutions obtained from P1dg~P2 applied to 
the /-plane. 

3.1. Discrete wave equation on the f -plane 



The P\bg-P2 spatial discretisation of the rotating shallow-water equations (see Cotter et al. 



(2009a) for a derivation) is 

^(/ U 5 } + WW X > = -c 2 <^,V^>, (2) 
±(4>',rf) = <V/V>, (3) 

where the velocity u 5 is in PIdgi t ne layer depth r] s = H{1 + r] s ) is in P2, for all test functions 
w s in PIdg an d (j> s in P2, and where c 2 = gH is the non-rotating wave propagation speed, g is 
the acceleration due to gravity, H is the mean layer depth and / is the Coriolis parameter. 

On the /-plane, / is a constant, and so we may take it outside the Coriolis integral. Applying 
the discrete Helmholtz decomposition to the velocity u s and the velocity test functions w s , i.e., 

u 5 = u 5 + V(j) 5 + V ± tp 5 + ii 5 , w s = w s + Va 5 + V ± f3 s + w 5 , 

equations become (after removing products of orthogonal quantities) 

(4) 
(5) 
(6) 
(7) 
(8) 

These solutions exhibit four types of orthogonal modes: geostrophic balance, inertia gravity waves, 
the physical inertial oscillation, and spurious inertial oscillations due to the presence of u. We 
shall now describe these modes one by one. 



A( Vc /,V/>-/ 


(Va\ V^> + c 2 (Va 5 ,Vr] 5 ) 


= o, 


d 

di 


(Va 5 ,W) + f(Va 5 ,V4> 5 ) 


= o, 






= o, 




A + 


= o, 




l t (w s ,u*) + f{w s ,(uY) 


= o, 
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3.2. Geostrophic balance 

For the continuous equations before discretisation, geostrophically balanced modes are obtained 



from non-zero steady solutions of the equations. As shown in Cotter et al. (2009a), in the PIdg-P 1 ^ 



discretisation solutions which satisfy the geostrophic balance relation are also exactly steady. To 
see this within the framework of this paper, assume a steady state, then equations (|4||8| become 



/ (Va s , W 5 ) + c 2 (Vc/, Vr] 5 ) 
/(Vc/,V0 5 ) 
-(Vc/,V0 5 ) 
f(w 5 ,u s ) 



0. 
0. 
0. 

o. 

0. 



(9) 
(10) 

(11) 
(12) 

(13) 



Equations (10) and (11) both imply that 0=0 since they are the usual continuous finite element 
discretisations of the Laplace equation which has no non-zero solutions because (f> s and a s are both 

ii 5 = 0. Equation Q is the 
and the Laplace operator can be inverted 



restricted to P2. Similarly equations (12) and (13) imply that u s 
discrete geostrophic balance relation between -^^and i] s , 
(since the finite element discretisation of the Poisson equation has a unique solution for solutions 
in P2) to obtain the pointwise geostrophic balance relation 



as noted in Cotter et al. (2009a). This means that P\dq-P2 has an excellent representation of 



,-.<5 



geostrophic balance. 

3.3. Inertia gravity waves 

The physical wave variables ip 5 and r/ 5 are uncoupled to the mean velocity component u 
and the spurious velocity component ii s . To obtain the discrete inertia gravity wave equation, the 
time derivative applied to equation Q gives 

(v^, W) - f± (W, v^> + 1/ (v^, V^> = 0. 



Substitution of equations ^ and (|6| then give 



+ f 2 )l t (* 5 ^) + l/(Va*,W)=0. 



(14) 



This is the usual continuous finite element discretisation of the inertia-gravity wave equation 



d 2 „o\ d 2 2 dr) 

Tj — C V 



dt 



+ r 



dt 



o. 



(15) 



Since only P2 functions are present, the solution i] s is third-order accurate, as opposed to the 
second-order accuracy expected with a first-order velocity discretisation. This higher-than-expected 
accuracy means that PI dg~P2 has a very accurate representation of inertia-gravity wave propa- 
gation. In particular, it should be expected that the phase velocity is more independent of mesh 
orientation than other second-order methods. The equivalent property for P0-P1 was noted in 



Roux et al. (2007), namely that the inertia-gravity dispersion relation was one order more accurate 
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than expected, namely second-order. The above proof extends this result to both arbitrary meshes, 
and to any finite element pair that satisfies the embedding properties above. 

A numerical verification of this third-order convergence is shown in Figure [Tj Care must be 
taken to obtain third-order convergence: if the initial conditions for the PIdg velocity are obtained 
by PIdg collocation, i.e. evaluating the analytic initial condition at the node points and using 
those values as nodal basis coefficients, then the truncation error in the initial condition for the 
velocity is second-order, and hence second-order accuracy is the most that can be expected after 
time-integrating the equations. However, a third-order accurate velocity initial condition can be 
obtained by first constructing a higher-order finite element approximation to the velocity field 
by collocation (we used a P2 approximation in the calculations in Figure [T]), and then applying 
the L 2 projection to obtain a PIdg velocity field. This results in third-order convergence of the 
free surface elevation over fixed time, since the free surface elevation equation is the P2 finite 
element approximation to the inertia-gravity wave equation, as shown above. To see that this 
procedure leads to a third-order accurate velocity field initial condition, first write the analytic 
initial condition for the velocity as 

u(x, 0) = V0 O + V ± ^o + wo- 

By standard approximation theory, the pth-order collocated finite element approximation to the 
initial condition satisfies u p = u(x, 0) + 0(Ax p+1 ). The P1 DG -P2 initial condition u s satisfies 

J v & -u s dV = J v s -u p dV 

for all PIdg test functions v s . After subsitution of the Helmholtz decomposition for u(x,0) and 
the discrete Helmholtz decomposition for u s , this becomes 

Va s -V(f> 6 dV = J Va s ■ V<p dV + 0{Ax p+1 ), 

jva s -V^ s dV = J Va 5 ■ Vip dV + 0(Ax p+1 ), 
u s = u + O{Ax p+1 ), 

and the potentials <f> s and ip s converge to 4>o and ipo as 0(Ax 3 ) following standard convergence 



theory for finite element discretisations of elliptic problems (see Brenner and Scott (1994), for ex- 



ample). Third-order convergence for the PIdg~P2 discretisation applied to inertia-gravity waves 



on the /-plane was demonstrated in Comblen et al. (2010) in which various partly-discontinuous 



finite element pairs were benchmarked against a high-order discontinuous Galerkin reference so- 
lution. Since the initial conditions were obtained by L 2 projection from the high-order solution, 
third-order convergence was observed. 

3.4- Physical inertial oscillation 

Since the integration is performed over spatially-independent functions, equation ^ may be 
written as 

w s ■ (f-u s + = 0, 

and since it must hold for all w s , we obtain 

which is the usual inertial oscillation equation which has spatially-independent solutions which 
rotate with frequency /. 
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x — x error (interpolated ICs) 

-- 0(dx~2 

x — x error (projected ICs) 

-- 0(dx~3) 



1.0e+00 



Figure 1: Plot showing convergence rates for the PIdq-P2 discretisation applied to the linear rotating shallow- water 
equations on an /-plane. The test problem is a single propagating sinusoidal wave in periodic boundary conditions, 
with the Li error in the free surface elevation computed after the wave has propagated all the way around the 
periodic domain. Second-order convergence for the free surface elevation is obtained when the initial conditions for 
velocity are obtained by collocation at node points; third-order convergence for the free surface elevation is obtained 
when the initial conditions are obtained by collocation with a quadratic P2 basis for velocity and reduced to the 
PIdg space by L2 projection. 



3.5. Spurious inertial oscillations 

Equation ([§]) describes the dynamics of the spurious velocity component it 5 . If 11 s 
velocity (i.e. is orthogonal to Va 5 and V^a* 5 ), then so is (u S )~ L and so equation 
involve any projection and hence can be written as 



is a spurious 
does not 



dt 



u 5 + f(u 6 



0. 



these solutions also simply rotate with frequency / and hence must be interpreted as spurious 
inertial oscillations which do not propagate as waves. 

If we replace the velocity space PIdg with the restricted space H(P2), as described in section 
[5J then we obtain the finite element pair which we call H(P2)-P2. we still have equations (|4]{7]) 
but without the spurious inertial oscillations in equation ([8]), hence the H(P2)-P2 discretisation 
has no spurious modes. 

3.6. Discrete dispersion relation for inertia- gravity waves 

In this section, we compute the discrete dispersion relation for the P\dg-P2 discretisation 
applied to the rotating shallow-water equations on the /-plane for the special case of a structured 
mesh in a regular hexagonal domain with edge length L centred on the origin, with periodic 
boundary conditions for opposing faces, tiled with equilateral triangles with edge lengths Ax = 
L/N for some positive integer N, and use this to define a continuous P2 finite element mesh. 



The discrete dispersion relation is developed by searching for time-harmonic solutions of (14). 
Assuming such a time-harmonic solution r/ s oc e tult , equation (14) becomes 



(-u 2 + f) (a 5 , V 5 ) + c 2 < Va 5 , V V 5 ) = 0. 



(16) 
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If r] S is an eigensolution of equation (16), then so is T z i] S (x) = r) S (x — z) for any z in the set V 



of translations that map vertices in the mesh to other vertices. Hence, eigenfunctions of equation 



(16) are all eigenfunctions of T z , i.e. they take the form 



Ikz 



$, Ax + z = x, VzgV, 



(17) 



where Q z is the translation of the hexagon formed from the six equilateral triangles surrounding 
the vertex at the origin by z, i7 <5 (^) is defined on the reference hexagon fl e with edge length 1 and 
centred at the origin, £ is the local coordinate in Q e , and k G IR 2 is the wave vector satisfying 
k ■ z = 2nl with I an integer. The wave vector k is contained in the first Brillouin zone of the 
periodic hexagonal domain which is bounded by the lines 



k ■ (cos(0 n ),sin(0 n )) J 



7T, 



9,, 



n + - I 7r/3, for n 



1,2,. ..,6. 



For more details of functions on periodic lattices, see (Kossevich, 2005), for example. 



Let us now fix an arbitrary wave vector k satisfying the conditions above. We note that the 



integral in equation (16) can be performed by integrating over all hexagons VL Z and dividing by 
three (since each equilateral triangle is covered by three hexagons). Given a test function a s , 



equation (16) (multiplied by three) becomes 



= [ {-u 2 + f 2 )a 5 (x)ri 5 (x) + c 2 Va s (x)-Vri s (x)dV(x) 

= J2 I ( Ax2 i-^ 2 + f 2 ) a5 U Ax + + V£c/(£ Ax + z) ■ e ik ' z d V{£), 

zev 

= / Ax 2 (-u 2 + f 2 )a s U)fj 5 ^) + c 2 Va 5 ^)-Vfj s U)dV^), 



where a s is defined on Q e with 

ze¥ 

We have now written the dispersion relation in such a way that all the computations can be done 
over one single reference hexagon Q e . The boundary conditions for i) s on the reference hexagon 
can be computed from the condition that rj 5 is continuous at the boundaries, meaning that on each 
edge of the hexagon fl e , denoted dfl e , n (with 1 < n < 6), 

where A£ is the vector from dQ e , n to the opposing face. Figure [2] illustrates the consequences of 
this for the basis coefficients of fj 5 when a nodal basi^j] is used. 



X A nodal basis is a basis in which each basis function has unit value at one of the node points, e.g. the vertices 
and edge midpoints in the case of the continuous quadratic mesh, and vanishes on all other node points. 
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Figure 2: Diagram showing the reference domain J7 e which is used to perform the numerical dispersion relation 
calculations. After considering the boundary conditions for fj s which are the consequence of requiring that r) S is 
continuous, there are four degrees of freedom for fj , which we denote {ij n }n=o- Each node in the diagram is labelled 
with a number n, indicating that i) s — fj n e lk '^ Ax at that node. 

We can similarly use continuity of a 5 to obtain boundary conditions for a <5 (^) on dfl e . On the 
boundary dfl ejn , 

= £V ({As - *)e**, 

= J2 a5 ^ + A €) Ax ~ ( z + ^xA$))e lk - z , 

zev 

= - A & Ax ~ z)e tHz - AxA t\ 

zev 

= e~ iAxk ^a\i + A£). 

This means that a s has boundary conditions which are the complex conjugate of the boundary 
conditions for fj S . 

We adopt a nodal basis for functions inside Vt e . There are 19 P2 nodes on VL e (see Figure [2J, 
and so we write 

19 

n=l 

where N n (g), (n = 1, . . . , 19), are the nodal basis functions for P2 functions inside fl e , and r\ n 
(n = 1, . . . , 19) are the nodal basis coefficients. The boundary conditions for ff described above 
can be expressed via a matrix S (which is a function of kAx due to the dependence of the boundary 
conditions for fj and a on fc), so that 

i) = Si), a = S*a, 
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where fj and a. are the vectors of the basis coefficients of rf and a s respectively, and fj and a. are 
the corresponding vectors of the independent degrees of freedom. 
After substituting, the wave equation becomes 



= (-oj 2 + f)(a 5 , V s ) + g(Va s ,Vrj s ) = 
= Ax^ T S^(-u 2 + f)M e + g ^jSf,, 



where M e is the local mass matrix 



and L e is the Laplacian matrix 

L e ,ij = [ VN^-VN^dV^), 

and f indicates the Hermitian conjugate of a matrix. Since ex. is arbitrary, we seek non-trivial 
solutions of 

St (Ax 2 (-co 2 + f) M e + gL e ) Si) = 0, 
and we obtain the dispersion relation 

\S^(Ax 2 (-u 2 + f 2 )M e + gL e )S\=0, (18) 

which must be solved for u given k (the k dependence is in S as described above). This equation 
is the determinant of a 4 x 4 matrix with entries that are linear in A = Ax 2 (u 2 — f 2 ), so it is 
quartic polynomial in A. 



After lengthy calculation using SymPy (SymPy Development Team, 2009), the following ma- 
trices are obtained: 



Sft M P S — { g T ^jj , Sft L e S — ( pT p 



where 

.4 



15 ' 



\/3 ^^cos(-ifc+|iV5) 

p _ ( reV3cos(ife+^V5) ^%/3cos(ifc) 

-iv^cos(iZ^3) -i^3cos(|fe-|iv^) 



c 

D 



15 ' 



^V3cos(|jfc+|l%/3) -^v^cos(fc)-^v^cos(|fc+|/ v ^)-^v / 3cos(-iA:+^v / 3) + |)V / 3 

8^ -| v / 3cos(-ife+iiV3) 
}-\/3cos(-ifc+ij%/3) 8^ 

-|%/3cos(|fc+|l-\/3) -fv^cosfAfc) 



, 3 y „™^ 4 „ , 4 „ v „y 3 v«-— ^ 2 'V nr|r | 

" " -fv^cos(ifc) _|V3cos(ife+i/^) ' ' 



F 



8^ -|^3cos(-|fe+j/v / 3) 
-|^cos(-|fe+li^) §v^cos(-|fc+|« v ^) + §v / 3cos(fc)+|v / 3cos(|A:+|«v / 3)+6v / 3 ' 
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having written k = (k, I). 

The resulting quartic equation for A = Ax 2 (u 2 — f 2 ) obtained from evaluating the deter- 
minant (18) is a very complicated expression that would take up several pages. Hence, solu- 
tions to the dispersion relation equation (18) were obtained by numerically evaluating the matrix 
(S^M e S)~ 1 S^L e S for various values of k, and using the Scientific Python linalg.eig routine, 
which were then sorted in numerical order. Since the equation for A = Ax 2 (u 2 — f 2 ) is quartic, 
this leads to four branches of the dispersion relation (this is typical for P2 schemes in two dimen- 
sions), which correspond to the fundamental exp(ifc • x) modes with k inside the first Brillouin 
zone, plus higher wave number solutions obtained from the second, third and fourth Brillouin 
zones which have the same translation property at the triangle vertices but result in different val- 
ues at the edge centres. The plots of the four branches are given in Figure [3] It is immediately 
visible that the lowest eigenvalues are very isotropic, as might be expected from the fact that 



the dispersion relation is in fact third-order rather than second-order, as described in section 3.3 



This means that resolved gravity waves of a particular wave number have a propagation speed 
which is largely independent of the direction of alignment of the mesh (this is a property which 
is considered important and was one of the contributing factors towards designing the hexagonal 
C-grid as an alternative to the triangular C-grid). It can also be seen that the dispersion relation 
is monotonically-increasing with |fe| with some small jumps when moving between branches (see 
Cotter et al. (2009b) for the equivalent one-dimensional plot); there are no spurious inertia-gravity 
modes. 



4. Discrete wave propagation on the /3-plane 

In this section, we consider the quasi-geostrophic scaling on the /3-plane, following the approach 
of Roux and Pouliot (2008); Thuburn (2008) in which the quasi-geostrophic approximation is 
applied to the spatially-discretised equations. 

In the /3-plane case, f = fo + (3y, and after substitution of the orthogonal decomposition for 
the solution variables and test functions into equations pj|3l we obtain 



dt 



(Vc/,V0 5 ) -/ (Vc/,V^) -(/3yVa s ,u s + (u 5 ) ± + W 5 + V V) + c 2 ( Vc/, Vr/ 5 ) = 



dt 



dt 
d_ 

dt 



(VoAW 5 ) + f (Va 5 ,V<p s ) + (PyVa 5 ,-u 5 - u s + V V + W> 5 ) = 



dt 



(c/,r/)-(Vc/,V0 5 ) = 0, 



<^n 5 ) + /o<^,(^) ± > + <^/3l/,(^) ± + (^) ± + VV-V^> = 0. 



At leading order in Rossby number in the quasi-geostrophic scaling, we obtain the geostrophic 
balance: 



-/o<V^V^} + c 2 (Vo/,V^) = 0, 

/o(Va*,V$) = 0, 

-(Vc/,V^) = 0, 

fo(w s ,(u s )f) = 0, 
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Figure 3: Plots showing contours of Ai(w 2 — f 2 ) in the kAx plane for each of the four branches of the numerical 
dispersion relation for the PIdc4-P2 finite clement scheme applied to the linear rotating shallow water equations 
on the /-plane. The lowest branch is shown top-left, with contours of the exact dispersion relation superimposed 
using dashed lines. This lowest branch is very accurate, and the contours are very circular, meaning that the wave 
propagation is almost independent of the direction of mesh alignment. The other plots show the higher branches 
which represent the second, third and fourth Brillouin zones in the kAx plane mapped in to the first Brillouin zone. 
For example, one can cross from the lowest branch into the branch in the top-right branch by going through the 
hexagon which bounds the region, emerging from the opposite edge in the hexagon in the top- right plot, moving 
in the opposite direction. It can be seen that all four branches represent physical modes from different regions of 
physical fc-space which can be resolved on the grid. 
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which we have already analysed in Section 3.2, and so we know that it implies that 



0, <# = 0, v 



/ 



V 



At the next order we obtain 
d 



dt 



Va 5 , V</4) - f < Vc/, V<> - (PyVa 5 , V^> + gH (Va 5 , VrQ 



dt 



09/ 



/o (w s , (u%) + -V^ 5 > 
fo (w s , (u%) + (w s /3y, -W g ) 






0. 
0. 





(19) 

(20) 
(21) 

(22) 

(23) 
(24) 



lag/ 1 \~ f {/> rg/ 

Notice that the spurious velocity modes do not appear at this order in the physical mode equations 
(20 22), and that equation (24) states that the ageostrophic spurious velocity modes are slaved to 
the geostrophic streamfunction. Substituting equations (19) and (22) into (21) gives 

Po 



_d_ 
di 



The second term in equation (25) may be written as 

(PyVa 5 , V^) = (VOW) - 0a* {0, 1), V"VJ) = ~P U 



0. 



<9 



(25) 



and we obtain the usual continuous finite element approximation to the Rossby wave equation 
using P2 elements 



d 

It 



0. 



(26) 



Since P2 elements are used, the approximation to the Rossby wave equation is third-order accurate, 
rather than the second-order accuracy one would expect with PIdg for velocity. The equivalent 
property for P0-P1 was shown in Roux and Pouliot (2008 ), namely that the Rossby wave dispersion 
relation is second-order. The above proof extends this result to arbitrary meshes and to any finite 
element pair which satisfies the embedding properties. 

We again expect that the phase velocity is more independent of mesh orientation than other 
second-order methods. Since the streamfunction ip s and the height variable r/ 5 are both from the 
P2 space and hence have the same numbers of degrees of freedom, there are exactly twice as 
many inertia-gravity wave modes as Rossby wave modes. We also note that if the reduced space 
H(P2)-P2 is used instead of P1dg-P2 we obtain the same equations but with vanishing spurious 
inertial modes. 

4-1. Discrete dispersion relation for Rossby waves 

Starting from equation (26), and following the method described above for the obtaining the 
inertia-gravity wave dispersion relation on the equilateral grid, we obtain the numerical dispersion 
relation 



itjj 



Ax 



S 



0. 



(27) 
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where D e is the local derivative matrix 



D 



and where / is the unit vector pointing in the direction of increasing / on the /3-plane. We shall 
investigate the variation in the dispersion relation with the alignment of the triangular grid, and 
hence it is convenient to write 



D e ,ij - h D l,ij + h D^u, 



where 



Jn 



After further algebraic manipulation with SymPy, we obtain 



f m^&dvit). 



SflDlS = 



Pi 



Qi 
Ri 



S^L P S = 



P2 

2 



Q2 
R2 



where 
Pi 

Qi 
Ri 



Pi = 



Ro 



-|v / 3sin(-ife+iiV3) 
-|v / 3sin(-ifc+i«v / 3) 

§V^sin(ifc+ijA/3) |V3sin(|fc) 

§>/3sin(§fc) -±V3sin(lk-llV3) + joV3sin(\k+llV3) 

-^v / 3sin(-|fc+|/ v / 3)- 1 i5V / 3sin(|fc+|/^3) 
-^V3sin(-ifc+|i v / 3)- 1 ^v / 3sm(|fe+|/v / 3) - g V3sin(k)- i V3sin( \k+\l^) - i a/3 sin ( Ifc-ljv^) 

| s in(-ife+i/v^) 
f sin(-|fe+|iV3) 

f sin( Jfc+IZv^) 
-isin(iiV3) - 1 ijsin(-|fe+iiv / 3) + ^sin(|fe+^v / 3) 





and 



'To sin (! fc - 



-^V3) + ^sin(-|fe+iZv / 3) 



fe+^%/3) + ^ sm(-\k+\lV2>) -yo sin(ifc+i/%/3)-^ sin(-ifc+§J\/3)- 



|jsin(|fc-|iV3) 



The eigenvalues can then be obtained using the method used for the inertia-gravity waves z. e. by 
finding the eigenvalues of the matrix for various kAx and plotting contours in k space. There is 
an extra difficulty in the Rossby case, because the numerical algorithm for obtaining eigenvalues 
of the 4x4 matrix does not preserve the order of the branches when kAx is varied. It is not 
possible to distinguish the branches by sorting the eigenvalues in numerical order for each k because 
the branches have values which cross. However, the branches can be distinguished by examining 
the corresponding eigenvectors. If we interpolate the continuous Fourier modes to the reference 
hexagon, we obtain four types of solution (after normalisation) for i), namely 



1 

4 
1 

4 



1 



1 

4 

^1 
4 



where the fundamental modes take the form of the vector on the left, and higher modes arise from 
the other three vectors. Hence, we identified the various branches by inspecting the eigenvectors 
and associating them with the branch which has the same sign pattern as the vectors above. 
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Figure 4: Contour plots showing w x 10 6 obtained from the solutions of equation (27 1, with parameters /o = 
1.0 x 10~ 4 , P = 1.0 X 10~ 12 , An; = 1.0 x 10 5 and c 2 = 1.0 x 10 5 (these parameters are the same as those used in 



Thuburn (2008)). / increases in the y-direction relative to the mesh. The lowest branch of the dispersion relation 



is shown top-left. The other branches are aliased higher values of kAx. 



Figure |4j shows contour plots of the frequency u for the case / = (0,1), with parameter 
values taken from Thuburn (2008). Exactly as the /-plane case, we obtain four roots for u which 
correspond to the fundamental modes (i.e. the modes that are possible to represent on a PI mesh) 
and the higher modes which arise from the extra accuracy on a P2 mesh. All the modes correspond 
to physical values after correct interpretation through the Brillouin zones as for the inertia-gravity 
wave case. A comparison with the exact dispersion relation for Rossby waves is given in Figure |5j 
a very close match is observed. Figure [6] shows contour plots for the same parameter values but 
with / = (1,0). Figure [7] shows the corresponding comparison with the exact dispersion relation; 
a close match is again observed. This shows that the P\dq-P2 discretisation has Rossby waves 
whose speed is almost independent of the mesh orientation. 



5. Summary and outlook 

In this paper we analysed the P\dq-P2 finite element pair applied to the rotating shallow- 
water equations, by means of a discrete Helmholtz decomposition which exists because of the 
embedding properties of PIog-P^i namely gradients and skew gradients of P2 map into PIdg- 
The discrete Helmholtz decomposition has some extra components, which we refer to as spurious 
velocity components, and which can be projected out, resulting in a discretisation that we referred 
to as H(P2)-P2. This decomposition was then used to show that in the /-plane, all steady 
states are geostrophically balanced (and vice versa). Furthermore, a discrete inertia-gravity wave 
equation can be derived which is the same as the P2 continuous finite element method applied 
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Figure 5: Comparison between the lowest branch of the discrete dispersion relation (left) and the exact dispersion 
relation (right). / increases in the y-direction relative to the mesh. 




2 2 4 



Figure 6: Contour plots showing w x 10 6 obtained from the solutions of equation (27 1, with parameters fo = 
1.0 x 10 -4 , P = 1.0 X 10~ 12 , Ax = 1.0 x 10 5 and c 2 = 1.0 x 10 5 (these parameters are the same as those used in 
Thuburn (2008)). / increases in the ^-direction relative to the mesh. The lowest branch of the dispersion relation 
is shown top-left. The other branches are aliased higher values of kAx. 
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Figure 7: Comparison between the lowest branch of the discrete dispersion relation (left) and the exact dispersion 
relation (right). / increases in the y-direction relative to the mesh. 



to the inertia-gravity wave equation, and hence the inertia-gravity wave solutions are third-order 
accurate. This should mean that the PIdg-P2 method should give very stable and accurate 
solutions of the linear geostrophic adjustment problem. We also showed that the spurious velocity 
components are uncoupled from the geostrophic balance or inertia-gravity waves, and they just 
undergo spurious inertial oscillations which do not propagate. When the H(P2)-P2 method is 
used, we obtain identical equations but without the spurious inertial oscillations. The H(P2)-P2 
method may be thought of as an implementation of the P2 finite element version of the Z-grid, in 
which vorticity, streamfunction and layer thickness are all collocated. Hence, the PIdg-P2 method 
may be thought of as a way to embed the finite element Z-grid into a method which avoids the 
need to solve elliptic problems for the potential and streamfunction, at the cost of adding spurious 
inertial oscillations. 

We then followed the methodology of Roux and Pouliot (2008); Thuburn (2008) to analyse the 
Rossby wave equation obtained from the the PI dg-P2 discretisation of the shallow- water equations 
on the /3-plane in the quasi-geostrophic limit. It was shown that the spurious velocity components 
do not couple in to the Rossby wave dynamics, in fact the geostrophic spurious components vanish 
and the ageostrophic components are slaved to the geostrophic streamfunction. It was shown 
that the quasi-geostrophic limit leads to a discrete Rossby wave equation which is identical to the 
continuous P2 finite element discretisation applied to the continuous Rossby wave equation, and 
hence the PIdg-P^ Rossby waves are third-order accurate. We expect that this means that the 
P\bq-P2 dispersion relation is much more independent of the direction of mesh alignment than 
other methods with linear velocity (such as the lowest-order Ravier-Thomas element which is the 
finite element version of the C-grid finite difference method). One seemingly negative aspect of 
using continuous finite element methods for pressure is that the mass matrix is not diagonal, so a 
linear system must be solved even when explicit timestepping is used. On the one hand, solving this 
linear system iteratively is extremely cheap since the condition number is independent of resolution 



and hence the number of iterations required stays constant under mesh refinement (Gresho and 



Sani, 2000). On the other hand, one can approximate the mass matrix M by a "lumped" diagonal 
mass matrix Ml with (Ml) a = ^2jMij. It was shown in (Le Roux et al. 2008) that lumping 



the mass has minimal effect on the dispersion relations so we would expect similar properties. In 
particular note that mass-lumping only effects the time-derivative terms so geostrophic states will 
remain steady. 

It seems almost inevitable (because of the difficulty in balancing the number of velocity and 
pressure degrees of freedom) that any numerical discretisation that is not based on quadrilateral 
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meshes will result in some form of spurious modes. From the results of this paper it appears that 
the PIdg-P'Z method puts the spurious modes into the least harmful place: it has no spurious 
pressure modes which would quickly pollute the solution and result in sub-optimal numerical 
convergence, it has no spurious Rossby modes which could modify the transfer of energy from 
barotropic to baroclinic modes in the presence of baroclinic instability, but it does have spurious 
inertial oscillations which do not propagate, and which can be filtered out using the H(P2)-P2 
projection. Whether or not these modes cause problems depends on how they are coupled to the 
physical modes through nonlinear advection, and this needs to be studied in careful benchmarks 
before recommending the P1dq-P2 method for use in NWP. If the modes are not harmful then the 
other properties discussed here (super-accurate wave propagation and representation of geostrophic 
balance on arbitrary unstructured meshes) mean that P1 DG -P2 should be an ideal choice for NWP 
models using adaptive mesh refinement. Here the projection filter will prove very useful, since the 
spurious modes can easily be extracted and measured, and modified advection schemes can be 
proposed which apply the projection before the wave step in semi-implicit splitting methods. 

Acknowledgements. This paper began after interesting discussions on spurious modes with Andrew 
Staniforth, John Thuburn and Nigel Wood. 
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